/* This function calculates two-dimensional interpolation. */

# include "scene.h"

extern short cmtrx_1dim_interpol(double *x, unsigned short nxpt, double *xipl, unsigned short nxiplpt,\
 unsigned short norder, struct cmtrx *cm, double **tempmtrx);


/* If the array **q is: q[nx][ny], the coefficient matrices should also be passed in the order
   of cmx, cmy. If this order is not maintained, the results will be errorneous.

   The icell, jcell speficied the cell-numbers at which the interpolation is to be done.
*/


extern short interpol_2dim(struct var_state *var, struct cmtrx *cmx, struct cmtrx *cmy, short icellx, short icelly, short icellz)
   {
    double qx, qy;
    short istartx, iendx, istarty, iendy, i, j, k, indxi, indxj, ireq, jreq, i2, j2, k2, dirz;

    ireq = icellx+(short)(2.01*var->dxreq);
    jreq = icelly+(short)(2.01*var->dyreq);

    qy = 0;
    istarty = cmy->istart[jreq-cmy->respadi];
    iendy = cmy->iend[jreq-cmy->respadi];

    k = icellz;
    dirz = 6 - var->dirx - var->diry; //dirx + diry + dirz = 1+2+3 = 6;
    for(indxj = istarty; indxj <= iendy; indxj++)
        {
         j = indxj+cmy->srcpadi-(int)(2.01*var->dydef);

         qx = 0;
         istartx = cmx->istart[ireq-cmx->respadi];
         iendx = cmx->iend[ireq-cmx->respadi];

         for(indxi = istartx; indxi <= iendx; indxi++)
            {
             i = indxi+cmx->srcpadi-(int)(2.01*var->dxdef);

             i2 = i*((4-var->dirx)/3) + j*((4-var->diry)/3) + k*((4-dirz)/3);
             j2 = i*(((var->dirx+1)%2)) + j*(((var->diry+1)%2)) + k*(((dirz+1)%2));
             k2 = i*(var->dirx/3) + j*(var->diry/3) + k*(dirz/3);

             qx += cmx->mtrx[ireq-cmx->respadi][indxi-istartx]*var->q[i2][j2][k2];
            }
          qy += cmy->mtrx[jreq-cmy->respadi][indxj-istarty]*qx;
         }

    var->qinterp = qy;
    return 1;
   }
